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ABSTRACT 

When an accretion disk surrounds a binary rotating in the same sense, the binary exerts 
strong torques on the gas. Analytic work in the ID approximation indicated that these torques 
sharply diminish or even eliminate accretion from the disk onto the binary. However, recent 
2D and 3D simulational work has shown at most modest diminution. We present new MHD 
simulations demonstrating that for binaries with mass ratios of 1 and 0.1 there is essentially no 
difference between the accretion rate at large radius in the disk and the accretion rate onto the 
binary. To resolve the discrepancy with earlier analytic estimates, we identify the small subset 
of gas trajectories traveling from the inner edge of the disk to the binary and show how the full 
accretion rate is concentrated onto them as a result of stream-disk shocks driven by the binary 
torques. 

Subject headings: accretion, accretion disks — binaries: general — MHD — methods: numerical 


1. INTRODUCTION 

Binary systems often form within a gaseous disk 
environment. Because the tidal torques exerted by 
the binary on the disk are repulsive when the disk 
and the binary orbit in the same sense, it has long 
been thought that these torques severely limit, or 
perhaps entirely prevent, accretion from the disk 


onto the binary Pringle (1991). As a result, a very 
low-density cavity forms within ^ 2 a of the bi¬ 
nary center of mass, where a is the binary’s semi¬ 
major axis. Despite this prediction, observations 
indicate accretion onto such binary systems at a 
level comparable to their single star counterparts. 
This is clearly detected for low-mass binary stars 


l nearby star-forming regions (e.g., White & Ghez 


2001). High velocity gas flows are observed bridg¬ 


ing gaps that are believed to be cleared out by pro¬ 
toplanetary companions (e.g., Casassns et al.|2013 


Rosenfeld et al.||2014 ). There are even detections 


of accretion onto planetary mass companions (e.g., 
Bowler et al. 2011 Zhou et al. 2014). A growing 
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number of dual AGN candidates have also been 
reported (e.g., Komossa et al.||2003 Rodriguez et 
al. 2006| Comerford et al. 2011); if systems like 


these become mutually bound, they could become 
accreting binaries. It is therefore important to un¬ 
derstand how gas is able to accrete despite these 
torques, both to be able to use electromagnetic sig¬ 
natures (Rodig ~et al.|2014 ) as diagnostics and also 
to learn how accretion influences the evolution of 
these systems. 

Much work has already been expended inves¬ 
tigating how the binary torque reshapes the sur¬ 


rounding circumbinary disk (e.g., Artymowicz & 


Lubow 

1994 

Crida et al.| {2006; MacFadyen & 

Milosav 

jevic 

2008 Shi et al. 2012] Noble et al. 


2012). However, how the torque regulates the gas 


accretion is a less well-developed subject. The ID 


analysis of Pringle (1991) has been extended, but 


adopting that paper’s assumption of zero accre- 


tion (e.g., 1 

Hilosavjevic & Phinney 

2005 Lodato 

et al. 2009 

Liu & Shapiro 

2010). More recently, 


there have been a number of 2D and 3D simula- 


tions of such systems (e.g., 

Artymowicz & Lubow 

1996 

Bryden et al. 1999; M 

ncFadyen & Milosavl- 

jevic 

2008 

Dotti et al. 2009 

Shi et al. 2012 Noble 

et al. 

2012 

D’Orazio et al. 2013 Farris et al. 2014). 
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Table 1 

Properties of Accretion Disk Simulations 


Label 

Type of Simulation 

Q 

Resolution a 

Radial Extent b 

(M(r in )) t c (VGMaEo) 

S2DE 

Hydrodynamics 

0.0 

256 x 64 

(0.8,40) 


S3DEQ 

MHD 

0.0 

512 x 400 x 96 

(0.8,40) 

0.0081 

S3DE 

MHD 

0.0 

512 x 400 x 384 

(0.8,40) 

0.0085 

B3DE 

MHD 

1.0 

512 x 400 x 384 

(0.8,40) 

0.011-0.014 

B3DEq 

MHD 

0.1 

480 x 400 x 384 

(1.02,40) 

0.013-0.017 


a Cell counts are listed as r x 0 x (f) for 3D MHD simulations and r x 0 for 2D hydrodynamic simulations. 

b ln units of a. The azimuthal extent is (0,27r) for all MHD simulations except S3DEQ, which is (0,7r/2). 

c Time 1 averaged accretion rate measured at the inner boundary. The accretion rate in the quarter-disk 
simulation S3DEQ is multiplied by 4. For binary runs, time averages at both early and late stages are 
presented and separated by a dash. 


These consistently find the cavity predicted ana¬ 
lytically, but also accretion at a rate comparable to 
or somewhat smaller than in the outer disk. Typ- 


single point-mass case; i.e. e 1 when q < 10 -3 . 
For values of q closer to unity, the results so far are 
somewhat mixed. Most simulations have found 


ically, the accretion takes place in narrow, high 

that e is reduced by a factor of a few (Bryden et 

velocity streams emanating from the edge of the 

al. 1999[|Lubow et al. 1999 

Shi et al. 12012] Noble 

cavity and spiraling inward toward the central bi- 

et al. 2012| D’Orazio et al. 

2013), but the depen- 

nary. Moreover, these studies have shown signifi- 

dence on q remains poorly defined. For example, 


cant consequences of the accretion for the develop¬ 
ment of the system. In protoplanetary disks, the 
amount of accretion through the gap controls the 


mass growth of a gap-clearing planet (Bryden et 
ah]|1999 Lnbow et al.||1999 D’Angelo et al.||2006 ) 


and the surface density within the gap (Zhuetal 
2011 Fung et ~aT7||2014 ). The angular momentum 
associated with the advecting gas also strongly in¬ 


fluence the orbital evolution of the binary (Shi et 
ah1|2012 Roedig et al.||2012 ). 


It is therefore crucial to quantify the net accre¬ 
tion fraction (e) of a prograde circumbinary disk, 
i.e. the ratio between the time-averaged net accre¬ 
tion rate M(q ^ 0) and the rate that would take 
place in the absence of binary torque M(q = 0): 


e ^M{q^ 0 ) 


(i) 


M(q = 0) 

For the extreme small mass-ratio case (we de¬ 
fine the mass-ratio q = M 2 /Mi < 1), previous 
studies using viscous hydrodynamics simulations 
have found that there is a threshold mass-ratio 
q ~ 10 -3 , such that for smaller q the accretion rate 
is nearly unchanged when compared to that for a 


D’Orazio et al. (2013) saw diminished accretion for 


q > 0.01 binaries, while Farris et al. (2014) saw no 


signs of suppression for 0.026 < q < 1. The situ¬ 
ation is further clouded by the fact that previous 
work (except for Shi et al. (2012) and Noble et 


al 


( |2012 )) has used the approximation of 2D hy¬ 
drodynamics, in which internal accretion stresses, 
both within the accretion disk itself and in the 
streams traversing the cavity, are described by the 
phenomenological a viscosity model, rather than 
the actual physical mechanism, correlated MHD 
turbulence. 

In this paper, we carry out 3D MHD simula¬ 
tions to investigate accretion in circumbinary disks 
around equal mass (q = 1) and unequal mass 
(q = 0.1) binaries. We assume the disk and bi¬ 
nary are coplanar, orbit in the same sense, and 
the binary orbit is a circle with fixed semi-major 
axis a. For a control experiment, we also simulate 
their single point-mass counterpart. Comparison 
between the single mass run and the circumbinary 
runs provides the answers to two key questions: 
(1) To what extent does the accretion rate of a bi¬ 
nary differ from that of a single object, i.e., what 
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is the value of e? (2) How does the accretion rate 
from circumbinary disks depend on the mass ra¬ 
tio q of the binary, i.e., how does e vary with g? 
As we will show, the answer to the first question 
is also the answer to the second: e ^ 1 for all 
q. This answer leads to a further question that 
we will also attempt to answer: what was missing 
from the early ID analytic studies that led them 
to conclude e^O? 

We organize this paper as follows: In §[2j we de¬ 
scribe the physical model and numerical setup of 
our single point-mass and circumbinary disk sim¬ 
ulations. In §[3j we present our simulation results. 
We analyze these results physically in §[4j Finally, 
we summarize our conclusions and discuss the im¬ 
plications of our findings in § [5] 


Numerical Simulations 


In this section, we discuss details about the nu¬ 
merical experiments we carried out in order to ex¬ 
plore the effects of the binary torques on the ac¬ 
cretion rate. We used the same code as IShi et ahl 
(2012), which is a modern version of the 3D, time- 
explicit Eulerian finite-differencing ZEUS code for 
MHD ( Stone fc Norman|1992a|b Hawley & Stone 


1995). 


2.1. Model Setup 


The model setup is very similar to the one de¬ 
scribed in Shi et al. (2012). Our account here is 


therefore very brief and emphasizes the few dif¬ 
ferences between our new simulations and the one 
presented in our earlier paper. We also summa¬ 
rize the main properties of all runs in this paper 
in Table [l] for reference. 

We construct our disks in an inertial frame with 
origin at the center of mass. We set the gravi¬ 
tational constant G and the central mass M to 
be unity, whether it is a single object or a bi¬ 
nary. The binary separation a = 1 is the sim¬ 
ulation lengthscal^] The time unit is therefore 
= y/a 3 /(GM)7 The density is normalized to 
the initial midplane value p 0 = 1 and the surface 
density unit is defined to be Eq = poa. Similar to 
Shi et al. (2012), we adopt a globally isothermal 
equation of state with a fixed sound speed, but 


1 In the q = 0 simulation, a has no physical interpretation 
other than the code-unit of length. 


we set c s = 0.1, twice that of the previous paper. 
A larger sound speed allows us to achieve better 
resolution and also perform longer duration simu¬ 
lations. Again, we assume the disk mass is consid¬ 
erably smaller than the mass of the central object 
and therefore neglect the self-gravity of the disk. 

As discussed in Shi et al. (2012, see section 2.2), 
the simulation grid needs to resolve three dif¬ 
ferent length scales: the disk scale height H , 
the maximum growth-rate wavelength of the MRI 
Amri = Stt/^/15va/^{R)^ and the spiral den¬ 
sity wavelength Ad ~ 27rc s /Ubi n , where va is the 
Alfven speed and Cl(R) is the disk rotational fre¬ 
quency. We constructed our grid in spherical co¬ 
ordinates following the same scheme as in 


Shi et 


al. (2012), which was originally proposed by 


Noble 


et al. (2010): logarithmic in the radial direction, 


uniformly spaced in azimuthal angle, and spaced 
according to a polynomial function in polar angle 
in order to concentrate cells near the orbital plane 
(Equation (l) of Shi et al. (2012), with £ = 0.825, 
0 C = 0.1 and n = 9). There were [AT r , 400,384] 
cells in (r, 0,0), covering a computational domain 
spanning [ri n ,r out ] radially, [0 c ,tt — 0 C \ meridion¬ 
ally and [0,27r] azimuthally, where r ou t = 40a, 
0 C = 0.1, ri n is the radius of the circular central 
cut-out (0.8a for q = 0 and q = l, increased to 
1.02a for q = 0.1 to confine the secondary inside 
the excision), and N r = 512 for q = 0 and q = l, 
diminished to 480 for q = 0.1 run as a result of the 
larger r- in . In terms of minimum cell size, our po¬ 
lar angle resolution is a factor ^1.4 better than 
in Shi et al. (2012). With a doubled disk scale- 
height (c s =0.1 vs. c s = 0.05), we therefore have 
about twice as many cells per scaleheight as in |Shi| 
et al. (2012). In recent years, standards have been 


developed for achieving resolution adequate to de¬ 
scribe the principal features of MRI-driven MHD 
turbulence: at least 10-20 cells per 2ttvaz,4>/Q, 
where VAz,(f> is the Alfven speed associated with 
the vertical (z) component of the magnetic field 
or azimuthal (</>) component ( Hawley et al.| 2Qll 


2013). We typically have ~ 20 cells per character¬ 


istic vertical wavelength and at least 10 cells per 
characteristic azimuthal wavelength. 

We choose strict outflow boundary conditions 
for the inner and outer radial surfaces and also at 
the edges of the polar axis cutout; all inward veloc¬ 
ities are set to zero on these boundaries. Periodic 
boundary conditions are used for the azimuthal 
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Fig. 1.— (a) Density distribution (color) at the end of the 2D hydro dynamical simulation S2DE, which is used to 
provide initial conditions for the 3D MHD quarter disk run S3DEQ. Note that the vertical scale is stretched relative 
to the horizontal. The black contour lines mark the initial magnetic field loops, (b) Density isosurfaces (color) of 
the bottom half disk of S3DEQ at t — 800^^. This data is used as the initial conditions of S3DE after combining 
four identical copies to make the whole 2tt disk, (c) Density isosurfaces (color) of the bottom half disk of S3DE at 
t = lOOOQ^. This data becomes the initial conditions for our circumbinary simulations. 
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direction. For the magnetic field in the radial 
and meridional directions, we set the transverse 
components of the field to be zero in the ghost 
zones. The components normal to the bound¬ 
aries are calculated by imposing the divergence- 
free constraint. 

2.2. Numerical Experiments 





Fig. 2.— Accretion history of the single mass quarter- 
circle simulation S3DEQ (black) and the full-circle 
S3DE MHD (red) disks. Top: Inner edge accretion 
rate as a function of time. Middle : Time-aver aged ac¬ 
cretion rate as a function of radius. Bottom : Mass 
interior to a given radius as a function of time. The 
accretion rates and enclosed mass for the quarter-disk 
are multiplied by 4 for comparison with whole disk 
values. Note the very similar behaviors between the 
quarter disk and the complete disk, in agreement with 
previous simulations: the time averaged accretion of 
MHD turbulent disks is nearly independent of the ex¬ 
tent of the azimuthal domain if it is more than tt/3 
( Hawley|2000 Papaloizou &; Nelson||2003 ). 


2.2.1. Single Mass Run: q = 0 

The single mass simulation both serves as a ref¬ 
erence point and provides the initial conditions for 
the circumbinary simulations discussed later. In 
order to speed up these lengthy simulations, we 
built a three-step ladder: 


Step 1. 2D Hydrodynamic Disk: We start 
from a 2D axisymmetric, inviscid, hydrodynamic 
simulation (S2DE) in the (r, 6) plane with the 
same numerical and physical parameters as de¬ 
scribed above except that the parameter £ in the 
#-grid definition is 0.92 and the number of cells 
is only lower 256 x 64 in r x 0. The goal of this 
2D simulation is to get rid of initial transients and 
to find a quasi-equilibrium disk solution. The ini¬ 


tial configuration follows that of Shi et al. (2012 


section 2.3): the disk stretches from r = 3a to 
r = 6a with constant midplane density po and 
is vertically hydrostatic. We evolve the disk for 
200011^, when the disk reaches a steady state 
(see Figure Qa)). 


Step 2. 3D MHD Quarter Disk: We then 
use the final dump from S2DE as the initial condi¬ 
tions for a 7r/2 wedge of a 3D MHD disk (S3DEQ). 
The initial disk for S3DEQ is therefore axisym¬ 
metric. Its initial density and angular velocity are 
interpolated from the results of S2DE data. Mo¬ 
tions in other directions are neglected as they are 
very small. The initial contours of magnetic vector 
potential A$ are a set of nested poloidal loops fol¬ 
lowing the contours of the density within the main 
body of the disk (the contour lines in Figure [lja)). 
The values of the A$ contours are A 0 (p — 0.1po)> 
where Ao is a constant determined by requiring 
the average plasma (3 = 100. The magnetic field 
B is computed by taking the curl of the vector 
potential. Run S3DEQ begins at t = 0 and lasts 
until t = 1330^1. A snapshot of the quarter-disk 
at t = 800D^ is shown in Figure [ljb). We find 
the quarter disk approaches steady accretion for 
r < 5a between t = 8000^ and 1300^7* (see 
Figure [2]). 


Step 3. 3D MHD full Disk: We then patch 
together four identical copies of the 7r/2 disk from 
S3DEQ at t = 8000^ to build the initial condi¬ 
tions for a full 27r disk (S3DE). This run continues 
from t = 8000^ through 134011^. The evolu¬ 
tion history (see Figure [ 2 ]), nearly the same as the 
quarter disk run, indicates that by ^ 800D7.* the 
full disk undergoes quasi-steady accretior0 We 


2 The time averaged accretion rate at different radii suggests 
our q = 0 disk is approaching inflow equilibrium for r < 4- 
5a (see the second panel in Figure [ 2 ]). Outside that radius, 
the accretion rate gradually shrinks, changing sign to out¬ 
flow at r ~ 7a. A small amount of mass outflow at large 
radius carries the angular momentum transported from in- 
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show a clipped isosurfaces plot of this fully turbu¬ 
lent disk at t = lOOO^in i n Figure [ljc). The disk 
at this moment serves as the initial conditions for 
the binary runs that are discussed next. 

2.2.2. Binary Runs: q 7 ^ 0 

We performed two binary simulations, one with 
mass ratio q = 1(B3DE) and one with q = 
0.1(B3DEq). For both runs, we use S3DE data 
at t = lOOOf^in as the initial conditions. Both 
also retain the domain size, numerical resolution, 
and physical parameters of the q = 0 run; the only 
change is to replace the point-mass at the center 
with either an equal-mass binary or a q = 0.1 bi¬ 
nary with the same total mass. In both binary 
cases, the orbit is constant and circular, and ro¬ 
tates counter-clockwise, prograde with respect to 
the disk. 

Everything about the q = 0.1 case is the same 
except for the extent of the radial grid. Because 
the secondary is ~ 0.9a from the center of mass, 
we increase the inner excision size from 0 . 8 a to 
fin = 1.02a. This change results in a truncation 
of the initial disk of the single mass run, but we 
keep all other aspects, such as the resolution and 
outer boundary of the domain, intact. 

3. Results 
3.1. q = 1 Run 

The initial disk adjusts to the new potential 
within 1-2 binary periods after the equal-mass 
binary replaces the point-mass. During this ini¬ 
tial transient phase, the binary clears out a cavity 
around itself in which the density is 10 -3 -10 -2 of 
what it was when there was a point-mass at the 
center. A small amount of mass initially found at 
r < 2a passes through the inner boundary (about 
2% of the total disk mass), but most of the disk 
mass within r < 2a is pushed outward by the bi¬ 
nary torque. After > 500^, the circumbinary 
disk gradually reaches a quasi-steady state. In its 
inner portion (r < 4a), the mass interior to r be¬ 
comes nearly constant in time (first panel in Fig¬ 
ure [ 3 ]) so that the surface density’s radial profile 
also changes only very slowly (Fig. [4]). We also 
note that the late-time average surface density of 


side. 




Fig. 3.— History of disk mass interior to a given 
radius for q — 1 (top) and q — 0.1 circumbinary disks. 


the outer disk follows the E oc r 2 
in 


Shi et al. (2012) quite well. 


scaling observed 



Fig. 4.— Surface densities for the initial condition 
(red dashed) and the time-averaged q — 0.1 (black) 
and q — 1 (green) cases. Solid curves are early-time, 
dashed curves are late-time averages. Solids are for 
the early time averages and dashed for the late time 
averages. For reference, a oc r~ 2 scaling is shown by 
the black dotted curve. 

After the first ^ 50D^, the accretion histo¬ 
ries show that the q = 1 disk has also reached a 
quasi-steady state with respect to this property 
(Fig. §. Clearly, the accretion is not hampered 
by the central binary (contrast the red dashed 
curve for the q = 0 disk with the green curve for 
q = 1). Although the strong binary torque clears 
out a low density cavity near the center, the cav- 


6 


















r _ __ 

- 



- t = 1050-1300, q = 0 

i 

- t = 1050-1250, q = 0.1 


- t = 1250-1500, q = 0.1 

Vn\ “ 

— t = 1050-1250, q = 1 


-- t = 1250-1400, q = 1 



1 10 
r / a 


Fig. 5.— Top : Accretion rates through the inner 
boundary of q = 0 (red dashed), q — 0.1 (black), and 
q — 1 (green) disks. Bottom : Time-averaged accretion 
rates for q = 0 (red dashed), q — 0.1 (black), and q — 1 
(green) as a function of radius. For the q ^ 0 cases, 
both early (solid) and late (dashed) time averages are 
presented. 


ity is not empty of gas. Streams of gas still flow 
through the potential maxima at the L2 and L3 
points (see Fig. 10) and maintain accretion at a 
rate comparable to the single mass case. The over¬ 
all time-averaged accretion rate of the q = 1 disk 
is in fact a bit greater than that of the q = 0 case: 
M(r = ri n ) ~ O.OH(GMa) 1 / 2 S 0 over t = 1050- 
1250(1^ (‘early-time’) and ~ 0.014(GMa) 1 / 2 So 
over 1250-14000^ (‘late-time’), compared to 
0.0085(GMa) 1 / 2 D 0 of the q = 0 disk (from 
t = 10000^* until the end of the simulation). 
These numbers translate to an accretion efficiency 
e ~ 1.3-1. 6 . This result is in approximate agree¬ 
ment with the viscous hydrodynamics simulations 
of D’Orazio et al. (2013), in which they found the 
accretion rate of an equal mass binary is ^ 0.93 
times the single-mass case (see their Table 3). It is 
in even better agreement with the results of Farris| 


et al. (2014), who found a ratio of 1.55. 


We show the radial dependence of the time- 
averaged accretion rate M(r ) in Figure [H] For 
both early (t = 1050-125011^) and late time av¬ 
erages (t = 1250-140011^), M(r) for both q — 1 


and q = 0.1 exceeds that of the q — 0 disk at 
all radii. The q = 1 disk achieves a fairly steady 
accretion within r < 4a, quite similar to a single¬ 
mass disk. At late times, the accretion rate at 
smaller radii increases by ^ 20%, and the zone of 
nearly constant accretion radius as a function of 
radius extends outward from r ~ 4a, characteris¬ 
tic of early times, tor^ 25a. 

The different levels of accretion observed in Fig¬ 
ure [5] are directly connected to variations in the 
internal stresses. As shown in Figure |6j although 
the Maxwell stress changes little with time, the 
Reynolds stress rises by a factor of ^ 2 after 
t = 12500^ at all radii outside r ~ 2 a. This 
increase in internal stress then drives a larger ac¬ 
cretion rate at late time. 

The increase in Reynolds stress occurs at the 
same time (^ 1250^4^) as the disk’s spiral den¬ 
sity waves change from a tightly-wrapped m — 2 
pattern to a single-armed wave (Fig. [ 7 ]). Although 
the spiral waves have little effect on the Maxwell 
stress, the vertically-integrated Reynolds stress is 
enhanced along the density crest of the m = 1 
wave. For instance, the region of large Reynolds 
stress between r = 4a and 8a at t = 1300£4kin 
(the yellow spiral arm winding from 9 o’clock to 6 
o’clock) is absent in the t = 1200^4^ snapshot. It 
appears that the single-armed spiral wave is more 
effective at conveying angular momentum outward 
than the two-armed wave. 

3.2. q = 0.1 Run 

The accretion flow in the q = 0.1 disk is quite 
similar to that seen in the q = 1 case. The 
disk reaches its quasi-steady state after the first 
^ 50£4^. The history of enclosed disk mass 
within a given radius in Figure [3] indicates a slowly 
evolving steady state in the inner part of the cir- 
cumbinary disk. Further evidence that an approx¬ 
imate state of inflow equilibrium has been reached 
for r < 4a is given by the time-averaged radial 
profiles of M(r) (Figure [ 5 ]). Also like the q = 1 
case, the accretion rate for q = 0.1 gradually in¬ 
creases over the course of the simulation, rising by 
- 30%. 

These two phases of accretion are also closely 
connected to the phase transition of the disk struc¬ 
ture for the q =5 0.1 run. In Figure |8j we show 
snapshots of disk properties before and after the 
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Fig. 6.— Time averaged stress to pressure ratios 
(top) and their associated angular momentum fluxes 
for the q — 1 disk. Different line styles denote dif¬ 
ferent time spans: t — 1050-1250 (solid), t — 1250- 
1400 (dashed), and t = 1050-1300 of the q — 0 run 
(dotted); different colors distinguish Reynolds (green) 
and Maxwell (red) stresses. There is an increase of 
Reynolds stress and its corresponding angular momen¬ 
tum flux at later times, owing to the m — 1 spiral 
waves. 


transition. It is clear that the disk evolves from 
a relatively compact two-armed spiral structure 
into a large scale single-armed structure. We no¬ 
tice that this transition happens faster than in the 
q = 1 run, possibly due to the asymmetry of the 
binary system itself. A common feature shared by 
these two phases is the gas stream attached to the 
secondary, which indicates stronger accretion onto 
the secondary than onto the primary. 

Like the q — 1 case, the accretion rate through 
the inner boundary in the q — 0.1 case (black solid 
curve in upper panel of Figure |5| appears to have 
two phases as well. At early times ( t = 1050- 
1250(2^), the time-averaged rate M(r = ri n ) ~ 
0.013(GMa) 1 / 2 So, while the accretion rate grad¬ 
ually increases so that the average at late times 
(t = 1250-1500) is ~ 0.017(GMa) 1 / 2 S 0 . In both 
phases, the time-averaged M(r = ?q n ) exceeds the 
q = 1 case by ^ 20%, and is ^ 1.5-2 times the 


single mass case [^] By comparison, in the 2D hy¬ 


drodynamic simulations of D’Orazio et al. (2013), 


the time-averaged accretion rate onto a q = 0.1 
binary was ~ 0.7x that of a single-mass system, 
while those of Farris et al. (2014) found a ratio 
~ 1.7. Just as for the q — 1 case, our results show 
a qualitatively similar but quantitatively some¬ 
what greater accretion efficiency than found by 


D’Orazio et al. (2013), and quite good quantita¬ 


tive agreement with Farris et al. (2014). 


3.3. Accretion rate fluctuations 


In both the q = 1 and q = 0.1 disks, the 
fractional amplitude of fluctuations in the mass 
accretion rate is a few times greater than for a 
point-mass (q = 0). As we will argue below, the 
larger fluctuations are driven by the binary itself. 
However, the accretion rate fluctuations in the bi¬ 
nary cases also appear to depend significantly on 
mass-ratio, both in amplitude and in frequency. 
The typical peak-to-trough amplitude contrast for 
q = 0.1 is ~ 0.01-0.015 Ho(GMa) 1 / 2 , about as 
large as the mean accretion rate. This is about 
twice the amplitude seen when q— 1. 

At early times in the q = 1 run, the accre¬ 
tion rate fluctuates periodically at a frequency of 


:e [9| 

’ by [IF 


This is the 
Orazio et al.l 


^ 1.5^bin, as shown in Figure 
stage, called the ‘transient state 
(2013), in which the disk is beginning to become 
elliptical and two streams of nearly equal strength 
run inward from the disk’s inner edge; |D’Orazio| 


et al. (2013) similarly found a periodic modula¬ 


tion with this frequency. Later in the simulation, 
the 1.5f2 bin peak in the power spectrum splits into 
two lower magnitude spikes, their frequencies cen¬ 
tered on the original peak frequency but separated 
by ~ O.lDbin* This split indicates a change in the 
disk structure from point-symmetric toward more 
eccentric shape. 


Unlike the q = 1 case, the q = 0.1 case can 
induce m = 1 disk asymmetry immediately. Con¬ 
sequently, right from the start we see strong peaks 
in the accretion rate power spectrum at 0.7Ubi n 
and 1.3Ubin* These modes are beats between 
the binary frequency Ubi n and the orbital fre¬ 
quency of the small density enhancement near the 


3 Here we have ignored the effects from the different ri n as 
they are negligibly small. The accretion rate of the equal 
mass binary hardly changes between r = 0.8a and r = 1.0a. 
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Fig. 7.— Four disk diagnostics before and after the spiral wave phase transition at ~ 12500^ in the q — 1 
simulation. Upper row : t — 120012^, lower row : t = 1300^^. From left to right, the variables shown are: surface 
density, density weighted radial velocity (normalized to the sound speed), vertically integrated Maxwell and Reynolds 
stresses in units of Soc^. The two white dots in the central cut-out region represent the binary members. White 
indicates off-scale high values. 
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Fig. 8.— The same diagnostics as in Figure |7| but for the q = 0.1 disk at t — 1100 and 12500^. In this case, the 
spiral wave transition took place at t ~ 1200fl^, slightly earlier than in the equal mass case. 
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Fig. 9.— Fourier decomposed power spectrum of the 
accretion through the inner boundary of a q — 1 binary 
(top panel) and a q = 0.1 binary (bottom panel). 


disk’s inner edge, ~ 0.3^bin- Another peak at ~ 
1.4£\in might be the first harmonic of the 0.7£\in 
beat. At later time, we still see the 0.7Dbin and 
1.4S7bin beat frequencies, but they shift slightly to¬ 
ward higher frequency as the disk gap expands by 
a small amount. Interestingly, we do not see the 
strong single peak at l.Oflbin observed by |D’Orazio 


et al. (2013) for this mass ratio. 


4. Analysis 

Having seen that e is actually slightly greater 
than unity, we now turn to an effort to under¬ 
stand this result. The first point to raise is that 
it is unlikely e > 1 can persist for a long period 
of time. If it were to do so, the inner region of 
the circumbinary disk (r > 2 a) would be drained 
of mass, inevitably leading to a reduction in the 
accretion rate onto the binary. Thus, the better 
way to think about the values of e seen in our sim¬ 
ulations, a few tens of percent greater than unity, 
is that the spiral waves excited in a circumbinary 
disk by the members of the binary create a suffi¬ 
cient enhancement of the Reynolds stress to raise 
the accretion rate per unit mass in the inner disk 
by a few tens of percent. By this means, an ac¬ 
cretion rate equal to that injected at large radius 



Fig. 10.— Left column: Time averaged midplane 
density (top), midplane (middle) and inner bound¬ 
ary (bottom) accretion rate pv r r 2 sin 0, both for the 
q — 0.1 binary over the last 5011^4 of the simulation. 
All figures in a frame comoving with the binary. Right 
column: same as left, but for q — 1 binary. Here neg¬ 
ative means inflow. The plus symbols in the midplane 
plots mark the L2 and L3 points. Summed separately, 
regions of inward and outward mass flux have compa¬ 
rable magnitude; their net, although smaller in mag¬ 
nitude, is consistently inward. 


can be sustained by a surface density somewhat 
smaller than required when the potential is due to 
a point-mass. Over longer times than we can fol¬ 
low with this kind of simulation, we expect that 
the surface density in the inner disk will decline 
to this level, leaving the disk in true inflow equi¬ 
librium. 

With that clarification, it is time to consider 
the question of why 2D and 3D simulations consis¬ 
tently see substantial accretion from circumbinary 
disks onto the central binary despite the contrary 
prediction made by ID studies. One clue to the 
answer comes from the structure of the accretion 
flow through the cavity: narrow streams. 

4.1. Stream Structure 

In the body of an accretion disk, the inflow 
speed is generically much slower than the orbital 
speed, ^ a(H/r) 2 v or \>, where a is the usual ratio of 
vertically-integrated stress to vertically-integrated 
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pressure, and H is the local scale height. On 
the other hand, this flow, although only ^ H 
thick in the vertical direction, takes place, on av¬ 
erage, around the entire circumference of the disk, 
through an area 2tt r wide. 

By contrast, the flow across the cavity (see 
Fig|To]) is restricted to very narrow streams. Along 
the central density maximum of the streams, they 
are typically ^ 2-3 H wide if measured sideways 
from the maximum to where the density drops by 
90%. Moreover, the density in the streams as they 
approach the inner boundary is ^ 3-10 times lower 
than the density in the disk body. Thus, in order 
to carry the same mass inflow, the inward veloc- 
ity in a stream must be ~ ( r/H)(p disk /p stleam ) 
larger than the typical disk inflow speed, or 
^ (pdisk/Pstream)cy(-F/r)u 0 rb* This condition can 
be easily achieved because the absence of stable 
closed orbits within r ^ 2a when q is not too small 
leads to a characteristic inward speed ^ 0.3u or b> 
whereas (pdisk/p s tream)a(#A) ~ 0.03-0.1 in the 
conditions of our simulation, and often smaller in 
real disks. 

Thus, one way of looking at the contrast be¬ 
tween the ID and the 2D/3D results is to observe 
that the ID picture, in which there is no inward 
motion within r ~ 2a, is correct at almost all, 
but not quite all, positions. We explore this idea 
further in the next subsection. 

4.2. Orbital properties of inward trajecto¬ 
ries 

The reason gas accretion occurs through such 
narrow channels is that only a small fraction of 
the gas possesses the proper initial conditions (in 
position-velocity phase space) for ‘infair trajecto¬ 
ries, i.e., trajectories that begin from the disk’s 
inner edge and reach the inner cutoff T[ n . Most 
of the gas near the disk’s inner edge that begins 
moving inward is turned around and flung back 
out to the disk by binary torques. This fate of the 
majority of the initially inflowing gas is what led 
the ID analysis to predict no net accretion. 

To test this explanation and identify exactly 
what the special conditions for successful inward 
flow are, we begin by solving a large number of 
restricted three-body problems further restricted 
to orbits entirely in the midplane. The binary 
for these trajectory integrations has a circular 


orbit and q — 1. The test-particles have ini¬ 
tial conditions evenly distributed within a phase- 
space volume defined by 1.5a < r < 2.5a, 0 < 
< j) < 2i r, —0 Ay/GM/a < v r < 0, and —0.4 < 
uJ — [flx(r) — Obin] < +0.4. Here is the angular 
frequency in a frame co-rotating with the binary. 
The initial specific angular momentum is therefore 
j = (u)' + fibin )r 2 . 

Our results are shown in the left panel of Fig¬ 
ure |Tl] in the form of a Parallel Coordinates plod 
When the initial radius is r = 2a (at which a cir¬ 
cular orbit requires j ~ lAh^/GMa when q = 1), 
particles whose initial j is > \.h\fGMa (a/Q^in > 
—0.63) cannot reach the inner boundary; particles 
with 1.3 < j/VGMa < 1.5 (-0.68 < ~ 

—0.63) can, but only if they also have v r < —(0.1- 
0.3 )y/GM/a\ particles with j < 1.3 y/GM a are 
able to travel to 7q n even with v r only slightly neg¬ 
ative. The behavior of both these classes of parti¬ 
cles can be understood by reference to the approxi¬ 
mate (i.e., ignoring the quadrupolar contribution) 
effective potential V e ^{r) ~ —GM/r + j 2 /{2r 2 ) 
at r = ri n . If v 2 < V e f r(rin), the particle radial 
kinetic energy is negligible, and infall can happen 
only when Wff (n n) < 0, i.e., j ~ 1.3 VGMa or less, 
or equivalently u/ is no more than — 0.68flbin* 
Alternatively, if v r is sufficiently negative, the par¬ 
ticle radial kinetic energy can be large enough to 
overcome a positive effective potential barrier. 

These conditions for infall are not easily met in 
a real disk because the mean j is close to the cir¬ 
cular orbit value, ~ 1.45\/GMa, and radial infall 
speeds < —(0.1-0.3 )y/GM/a are rare in the disk 
body. Consequently, only those few fluid elements 
with j well below the mean can contribute to the 
inflow, and they are found in a tightly-constrained 
portion of phase space. In the right panel of Fig¬ 
ure ED we show a selection of fluid elements taken 
from a snapshot of our simulation at t = 120811^. 
The samples are drawn from a ring of disk around 
r = 2a with a width of 0.5a in order to capture 
the inner edge of the disk. Comparing the orange 
regions (all cells) with the blue ones (those falling 
in), one can conclude that indeed only a tiny frac¬ 
tion of the disk proper falls to the binary, while 
the rest of its fluid elements are pushed back out 


4 Parallel Coordinates is a common method for multidimen¬ 
sional visualization. For our case, each particle is repre¬ 
sented by a line connecting multiple attributes shown as 
vertical coordinates in parallel. 
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Fig. 11.— Left: A parallel coordinate plot for test particles at r = 2a distributed uniformly in 0, v r (shown in units of 
(GM/a ) 1 / 2 , and uJ — Q — flbin- The last coordinate j = (u/ + flbin)^ 2 is the initial specific angular momentum in units 
of (GMa) 1 / 2 . The green lines show the initial locations of all infalling particles; note that all have j < 1.5(GMq) 1//2 , 
the specific angular momentum that supports a circular orbit in the region sampled. The red lines show particles 
with enough angular momentum to trace nearly-circular orbits, but also fall inward; all have especially large inward 
speed as well. The red particles concentrate at two opposite azimuth angles where the binary potential is relatively 
weaker than at other angles, making it easier for particles to fall in. 

Right: Similar to the left but now the lines represent individual cells drawn from the actual simulation of the q — 1 
disk at t — 1208. The orange lines show all cells within annulus having r ~ 2a; the blue lines show all the infalling 
cells; the green lines show infalling cells with initial radius exactly 2 a. Comparing the orange and blue regions, only 
a small fraction of the disk cells actually fall in. The green cells found here are consistent in their properties with 
the green test particles shown on the left: v r <0 for u/ < —0.68. 


even if they initially move in a small distance. The 
criterion governing which fluid elements are able 
to move inward is identical to that identified for 
the test-particles: j < 1.3a/ GMa. In contrast, 
most of the orange fluid elements have either too 
large a j (or a/), i.e. large f/ff, or too small an 
inflow velocity — v r , and are therefore flung out. 
Particularly low j fluid elements move inward so 
quickly there is too little time for binary torques to 
substantially raise their angular momentum; be¬ 
cause they can pass the effective potential barrier 
at ri n , they cross the inner boundary and are per¬ 
manently removed from the circumbinary disk. 


Because gas flows near the binary are close to 


ballistic (Lubow & Artymowicz 2000 Shi et al. 


2012), we can estimate the accretion rate by count¬ 


ing the number of fluid elements satisfying the in¬ 
flow criterion. The data shown in the right panel 
of Figure [TT] indicate that the blue fluid elements 


cover ~ 1/10 of the annulus between r ^ 1.9a and 
2.2 a. For a typical surface density ~ 0.2H 0 and an 
infall timescale ^ 27r/f2bin, the inferred inflow rate 
would be M ~ 0.012v / §MaS 0 , consistent with 
our simulation results. Thus, although 90% of in¬ 
ner disk fluid elements cannot go any significant 
distance inward, the angular momentum distribu¬ 
tion is broad enough that its low j tail suffices to 
carry the full accretion rate. 

Having found that the fluid elements able to 
accrete are defined by their low specific angular 
momentum, the next question to answer is how 
their angular momentum is reduced to that level. 
Ordinary MHD turbulence does not broaden the 
angular momentum distribution to this degree: as 
shown in Fig. [6j the Reynolds stress in the q = 0 
case is only 1% of the pressure, and the pres¬ 
sure is only ~ 1% of the orbital energy per unit 
mass. Instead, the answer appears to be a conse- 
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quence of the binary torques themselves. As pre¬ 
viously remarked, most of the mass in the streams 
that move inward from the circumbinary disk re¬ 
turns to the disk after its specific angular momen¬ 
tum is raised by those torques. With that addi¬ 
tional angular momentum, its azimuthal velocity 
is somewhat greater than the local orbital veloc¬ 
ity when it strikes the disk (j ~ 1.6(GMa) 1 / 2 
as opposed to j ~ 1.45 (GMa) 1 / 2 ). The work 
done by the torques also increases the streams’ en¬ 
ergy, giving them an outward radial speed ~ (0.3- 
0 ^(GM/a) 1 / 2 . 


log 10 (Z/Z 0 ) J ™ 5 _ 1>0 _™ 0 



Fig. 12.— Surface density (color contours) and mass- 
weighted vertically-averaged velocity in the orbital 
plane (arrows) at t— 1209.6. 


An example of such a rapidly-moving outward 
stream can be seen near r ~ 2.3, <j) ~ 1.5 in Fig- 
Part of the stream (r ~ 2.3-2.4, ^ 1- 
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ure 

1.5) is still heading outward. However, a part of 
that stream has already encountered the ridge of 
high density trending from larger radius to smaller 
between 0 ~ 2 and ~ 3. When it struck that 
ridge, a pair of shocks was formed, a forward shock 
propagating into the disk material and a reverse 
shock propagating into the material that had been 
moving outward. In the reverse shock, a portion 
of the stream mass loses a significant fraction of 
both its angular momentum and its orbital energy 
and heads back inward, seen in this snapshot along 
the line from (r ~ 2.3, </> 2) to (r ~ 1.8, 0 ~ 3). 


This shock deflection is the origin of the low angu¬ 
lar momentum material which can now travel all 
the way to the inner boundary. 


Because this mechanism depends almost en¬ 
tirely on 2D orbital mechanics, its nature should 
be only weakly dependent on parameters such as 
the disk’s aspect ratio H/r = c s /v or b, provided 
only that c s /v or b is small enough that the stream- 
disk interaction is supersonic. In our simulation, 
in which the sound speed is 0.1(GM/a) ly/2 , the 
Mach numbers of these shocks are ^ 3-5; in real 
disks, they could be considerably larger. Support 
for the view that the stream-disk interaction is at 
most weakly affected by the disk thickness also 
comes from the fact that 2D simulations employ¬ 
ing laminar hydrodynamics and an “a viscosity” 


to mock up the fluid’s internal shear stress (Far¬ 


ris et al. 2014) find a very similar continuity in 


the accretion rate. It is possible, however, that 
the dynamics of the stream-disk shocks may de¬ 
pend upon the equation of state of the shocked 
gas. In our simulation, we assumed that the gas 
is isothermal. If, instead, this sort of shock were 
to occur in a less lossy gas, the immediate post¬ 
shock temperature would be much greater than 
the disk gas temperature. The shocked gas could 
then swell vertically well beyond the thickness of 
the disk, permitting it to flow above and below the 
disk. Investigation of such effects is well beyond 
the scope of our effort here. 


4.3. One-armed Spiral Wave 

As we have already remarked, a surprising out¬ 
come of our simulations is that at late times 
a strong one-armed spiral wave propagates out¬ 
ward through the circumbinary disk, enhancing 
the Reynolds stress sufficiently to increase the ac¬ 
cretion rate by a few tens of percent. Its origin is 
worth further attention. 

The most likely cause of this m = 1 feature 
is the complement of the mechanisms studied in 
the previous subsection. There we focused on the 
properties of those special fluid elements able to 
travel all the way from the inner edge of the disk 
to the binary. Here we focus on all the others, the 
ones that may initially move inward, but then feel 
the torques exerted by the binary and are pro- 


5 Also available at http://www.astro.princeton.edu/ 
"jmshi/leakage.htm 
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pelled outward again until they strike the disk’s 
inner edge. To demonstrate how these streams 
excite spiral waves, we fit the wave form seen in 
our simulations (surface density plots in Figure [7] 
and[8| with a standard density wave pattern (e.g., 
Equation (36) of Rafikov |2002| ). The resulting pat¬ 
tern speed is ~ 0.17f^bin (— 0.21Dbin) for the q = 1 
(q = 0.1) case, which suggests wave excitation oc¬ 
curs at r 3.2a when q — 1 and r ~ 2.8a when 
q 0.1. As shown in Figure [7] and [8j the density 
is strongly enhanced near r ^ 3a. This connection 
is most clearly seen in an animation (Fig. [I3|). 

The animations further show that the point of 
impact, the azimuthal location of the spiral wave 
driving point, rotates around the disk’s inner edge 
at an angular frequency slightly lower than the 
binary frequency. 


Similar to tidally or mass-transfer induced spi¬ 
ral shocks that transmit angular momentum out¬ 
ward in circumstellar disks in a close binary sys¬ 
tem (e.g., Sawada et al.|[l9 86 ; Rozyczka & Spruit 
1993), the waves excited by the binary-driven 


streams can also provide long-range angular mo¬ 
mentum transport. Figure [l4| shows the net out¬ 
ward angular momentum flux (AMF) associated 
with the time averaged Reynolds stress before and 
after large scale m = 1 density waves are excited 
(see green curves in the bottom two panels). Com¬ 
pared to the early period (top row), the Reynolds 
stress contributes a sizable negative torque at disk 
radius ^ 3-4.5a (2.5-4a) at late time (bottom row) 
in the g = 1 (g = 0.1) disk, which could explain the 
slightly greater accretion rate observed in q ^ 0 
disks than in the q = 0 disk. 


These large scale one-armed spiral density fea¬ 


tures were not previously reported in Shi et al. 
(2012), D’Orazio et al. (2013), or |Farris et al. 


(2014k We believe they are due to the greater 
global isothermal sound speed (O.lflbin^) used 
here, twice the sound speed in simulation B3D |Shi| 
et al. (2012). The faster sound speed causes spiral 


waves to stretch farther radially, creating large- 


al. 1994 

). In D’Orazio et al. 

(2013| (as well as in 

(MacFadyen & Milosavljevic 

2008)), H/r = 0.1 at 


all radii, so that the isothermal sound speed drops 
oc r -1 / 2 . The wavelength becomes shorter as the 
wave moves outward, causing the wave to become 
more and more tightly wrapped as it travels to 
greater radius. Given this apparent sensitivity to 




Fig. 14.— Radial derivatives of time-averaged and 
shell-integrated angular momentum flux (AMF) as 
functions of radius for the q — 0.1 (left) and q — 1 
(right) cases before and after the phase transition, 
where red shows the differentiated Maxwell AMF and 
green represents the differentiated Reynolds AMF. 
Negative means outward transport of angular momen¬ 
tum. Time averages are long enough that the net 
change of local angular momentum are close to zero. 
The late time averages of the q — 1 (q = 0.1) cases 
show that the Reynolds stress dominates the angular 
momentum budget at r ~ 3-4.5 a (2.5-4a), a location 
coinciding with the stream impact regions shown in 
Figure [7] for the q— 1 case and Figure [8] for the q — 0.1 
case. 


the disk equation of state, it remains to be seen 
whether such strong spiral waves are generic or 
rare in real circumbinary disks. 


5. CONCLUSIONS 


In this paper, we carried out 3D global MHD 
simulations of circumbinary disks with binary 
mass ratios 0.1 and 1 and contrasted them with a 
disk orbiting a solitary point-mass. We found two 
major results: 


1. The time-averaged accretion rate from a cir¬ 
cumbinary disk with either q = 1 or 0.1 
is indistinguishable from that of a circum- 
solitary disk whose central mass is the same. 
In other words, essentially all the mass sup¬ 
ply given the disk at large radius ultimately 
leaves its inner edge and travels to the bi¬ 
nary. The similarity of the q = 1 and q = 0.1 
cases in this regard suggests that this re¬ 
sult depends at most weakly on binary mass- 
ratio. This result confirms, with physical in¬ 
ternal fluid stresses, the conclusion reached 
by Farris et al. (2014) on the basis of 2-d 
hydrodynamics and a phenomenological vis¬ 
cous stress. 


2. The key reason why initial ID analyses sug¬ 
gested that the accretion efficiency (the pa- 
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rameter we call e) is ~ 0 rather than ~ 1 
is that they omitted consideration of the 
small volume in orbital phase space from 
which trajectories can travel inward from 
the disk’s inner edge all the way to the bi¬ 
nary, avoiding the strong torques that, for 
most of phase space, push matter back out¬ 
ward. Only those fluid elements with spe¬ 
cific angular momentum ~ 15% less than 
the circular orbit value at the disk’s inner 
edge can cross the gap and reach the bi¬ 
nary. Ironically, fluid elements with exactly 
this property are created as a consequence of 
the binary torques themselves: these torques 
add enough angular momentum to other gas 
traveling through the gap to propel it back 
out to the disk; it shocks upon reaching the 
disk, and a portion of its mass is deflected 
onto accretion orbits. 


Our results have several observational implica¬ 
tions. If all the matter supplied to a circumbi- 
nary disk at large radius accretes onto the binary 
through narrow streams, those streams must shock 
when they strike the outer edges of accretion disks 
around the members of the binary. In the context 
of the formation of stellar binaries, these shocks 
may be the sites of the v = 1-0 5(1) H 2 vi¬ 
brational lines that can often be detected (Beck 


et al. 2012). Radial inflow streams at velocities 


approaching free-fall can also potentially explain 
both the observed low density cavities in transi¬ 
tional disk systems and the relatively normal stel¬ 


lar accretion rates they maintain (Rosenfeld et al. 


2014). In the context of supermassive black hole 


binaries, one can similarly expect strong shocks 
where the streams strike the outer edges of the 
individual disks. Hard X-rays from those shocks 
may be observable when the binary separation is 


small enough (Roedig et al. 2014 Farris et al. 


2015). If the accretion flow fed in at large radius is 


large enough, still greater luminosity can be gen¬ 
erated when the accreted matter approaches the 
two black holes’ IS CO regions. 
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Fig. 13.— Streams returning to the circumbinary disk exciting spiral waves. Click the figure to play a short 
animation of surface density in the q — 1 and 0.1 cases |^] Spiral waves can clearly be seen to begin when a stream 
pushed outward by the binary torques strikes the inner edge of the disk. 
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